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ABSTRACT 

For planets other than Earth, particularly exoplanets, interpretation of the 
composition and structure depends largely on comparing the mass and radius 
with the composition expected given their distance from the parent star. The 
composition implies a mass-radius relation which relies heavily on equations of 
state calculated from electronic structure theory and measured experimentally 
on Earth. We lay out a method for deriving and testing equations of state, and 
deduce mass-radius and mass-pressure relations for key, relevant materials whose 
equation of state is reasonably well established, and for differentiated Fe/rock. 
We find that variations in the equation of state, such as may arise when extrap- 
olating from low pressure data, can have significant effects on predicted mass- 
radius relations, and on planetary pressure profiles. The relations are compared 
with the observed masses and radii of planets and exoplanets, broadly supporting 
recent inferences about exoplanet structures. Kepler- 10b is apparently 'Earth- 
like,' likely with a proportionately larger core than Earth's, nominally 2/3 of the 
mass of the planet. CoRoT-7b is consistent with a rocky mantle over an Fe-based 
core which is likely to be proportionately smaller than Earth's. GJ 1214b lies be- 
tween the mass-radius curves for H 2 and CH 4 , suggesting an 'icy' composition 
with a relatively large core or a relatively large proportion of H 2 0. CoRoT-2b is 
less dense than the hydrogen relation, which could be explained by an anoma- 
lously high degree of heating or by higher than assumed atmospheric opacity. 
HAT-P-2b is slightly denser than the mass-radius relation for hydrogen, sug- 
gesting the presence of a significant amount of matter of higher atomic number. 
CoRoT-3b lies close to the hydrogen relation. The pressure at the center of 
Kepler- 10b is TPa. The central pressure in CoRoT-7b is probably close 

to 0.8 TPa, though may be up to 2 TPa. These pressures are accessible by pla- 
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nar shock and ramp loading experiments at large laser facilities. The center of 
HAT-P-2b is probably around 210 TPa, in the range of planned National Ignition 
Facility experiments, and that of CoRoT-3b around 1900 TPa. 

Subject headings: Equation of state - Planets and satellites: composition - Planets 
and satellites: interior 
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Introduction 



Planets outside the solar system have been detected since 1992 (IWolszan fc Fraillll992h 
from the Doppler shift of spectral features in emission from the star, which determines 



the o rbital period and places a constraint on the mass of the planet ( IMayor fc Queloz 



19951 ). Since 1999, the presence o 



across the face of the parent star (IHenry et al 



exoplanets has a lso been deduced from their transit 



20001 ) . The fraction of light blocked by the 



planet allows the radius o 



( ICharbonneau et al 



2000 



the planet to be deduced as a function of the radius of the star 



Carter et al 



20111 ). Subsequently, several hundred exoplanets 



have been detected at different distances from their stars, and the precision with which 
mass and radius have bee n deduced has i ncreased for some exoplanets to better than 10% 



in both mass and radius ( Schneider 



20111). 



20021) 



In 2001, light from an exoplanet was detected directly (ICharbonneau et all 
opening the window to studies of exoplanet structure through the composition and 
temperature of the surface or atmosphere. However, inferences about the composition and 
structure rely on the comparison of mass and radius with planets within the solar system. 
With the exception of Earth, which is currently the only planet for which seismic data 
exist, interpretations of the internal structure of the planets rely in turn on assumptions 
about the composition and temperature profiles through the planet. 

Theories of planetary formation can be investigated by comparing the structure of 
exoplanets with those within the solar system. Another motivation is to estimate the 
occurrence of Earth-like planets, in terms of mass and composition, and also those that 
might be habitable for life. 

Deductions about planetary structures, i.e. the composition profile, depend on 
the compressibility of the possible compositions thought to occur. The compressibility 
is needed over the range of pressures and temperatures occurring within each planet. 
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The compressibility is found from the derivative of the pressure-density relation at the 
appropriate temperature, which can be determined from the equation of state (EOS) 
for the particular composition of matter of interest. The development of EOS has been 



Mc Queen et al 



2001 



1970) 



Stacey fc Davis 



geophysics (e.g. 



2004 



Ahrens fc Gregson 



Stixrude et al 



1964 



impacts 


(for 


Ahrens 


1966 



Alfe et al 



2009), and inertial confinement fusion (fLjndl 



19981 ). There is a perception that experimental and theoretical methods for determining 



EOS are not available in regimes necessary to understand the internal structure of planets 



for pressures between 200 GPa and 10 TPa, i.e. fro m the limit o 



the onset of the T homas- Fermi-Dirac (TFD) regime (ISeager et al 



diamond anvil data to 



2007 



Grasset et al 



Some studies (e.g. 



Seager et al 



2009|). 



20071 ) have considered sub-TFD EOS with care, but it is 
common practice even when accurate theoretical calculations are available to represent the 
material with ad hoc functional forms, potentially leading to problems when extrapolating 
beyond the bounds of the constraining data and in switching over to TFD at high pressures. 

Although there is a definite need for more theoretical and experimental studies, 
appropriate theoretical techniques are just as accurate above 200 GPa as below, and are 
more than adequate to draw inferences about the structure of exoplanets; moreover, shock 
and ramp loading experiments can readily explore states up to ~5 TPa and are suitable to 
test and calibrate EOS. In this paper, we study the validity of electronic structure methods 
for predicting EOS up to the o(100) TPa pressures Q apparently occurring in exoplanets, 
and the capability of dynamic loading experiments to measure relevant states. We calculate 
mass-radius relations for several compositions of matter representing different classes of, 
and layers in, planets, discussing the validity of the EOS used. Finally, we compare the 
mass-radius relations with representative planets and exoplanets. 



l A pressure of 1 TPa is 10 million atmospheres. 
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2. Matter at high pressure 

The quasistatic structure of self-gravitating bodies depends on the scalar EOS, which 
can be described by an appropriate free energy expressed as a function of its natural 
variables, such as the Helmholtz free energy f(p,T), where p is mass density and T 
temperature. In principle, one could consider the composition as a set of explicit parameters 
in the EOS. In practice, planetary structures are thought to comprise layers in each of 
which a single composition, or range of compositions, dominates the EOS, such as Fe in the 
core of rocky planets. Therefore, we consider a separate EOS for each layer. 

As with dynamic loading situations, the pressure p = Of /Ov\t is the most directly 
useful aspect of the EOS for calculating and interpreting planetary structures. Almost 
invariably, a thermodynamically incomplete EOS is used: p(p, T) or p(p, e) where e is the 
specific internal energy. Planetary structures may be constrained to follow a specified 
temperature profile, in which p(p, T) is the more convenient form, or an isentrope, for which 
p(p,e) is convenient since p = —de/dv\ s where v = 1/p. 

In planets, brown dwarfs, and main sequence stars, the EOS is dominated by 
electrostatic forces and Pauli exclusion among the electrons and ions, rather than by 
strong-force interactions between the nuclei. In stars, the radiation pressure must be 
included, and the entropy is much higher, exploring a different region of the EOS. 

2.1. Theory 

In this section, we describe key theoretical methods and approximations used when 
predicting EOS, and comment on their applicability to states likely to occur in exoplanets. 

The relevant EOS of matter can be calculated using electronic structure theory. For 
a representative set of atoms defining the composition, thermodynamic potentials can 
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be calculated as a function of p (i.e. system volume) and T. Because the mass of the 
electrons is so much less than that of the nuclei, the state of the electrons can generally be 
considered with respect to the instantaneous positions of the nuclei - the Born-Oppenheimer 
approximation. Forces on the nuclei can also be calculated with respect to their positions 
and the distribution of electrons. The energy of the system can be thought of as comprising 
the ground state energy for stationary ions, thermal motion of the ions, and thermal 
excitation of the electrons out of their ground state. In principle, all contributions should 
be calculated self-consistently. However, for EOS, it is generally a good approximation to 
partition the free energy into the cold compression curve, thermal motion of the ions, and 
thermal excitation of the electrons, 



f(p,T) = f c (p) + f t (p,T) + f e (p,T) 



(1) 



sec 



1 



Swift et all 120011 ) . To describe the material state and the dynamics of most atoms, the 



quantum nature of the nuclei can be ignored, allowing their treatment as point-like masses 
exhibiting Newtonian dynamics. The exception is H, for which the quantum behavior of 
the proton can be represented as a correction to the potential experienced by the electrons. 

The motion of the nuclei can be calculated with respect to time. This technique, 



known a s first-principles molecu 



(QMD) (IMichielsen k De Raedt 



ar dy namics (FPMD) or quantum molecular dynamics 



19961 ). is particularly appropriate for fluids, unknown or 



ill-defined crystal structures, and multi-species compositions of matter where it is not clear 
a priori where particular species will be located. In addition to the EOS, FPMD has 
been used to predict transport properties, including the electric and the rmal conductivit ies 



( Recoules et al 



Clerouin 



2009 



Hamel et al 



2011 



Hoist et al 



2011 ) and viscosity f Alfe et a 



2000a; 



20021 ). which are important in understanding the formation and thermal profile 



of planets as well as the generation of their magnetic fields. FPMD has also been used in 
studies of the free energy profiles of mixtures to determine possible miscibility gaps which 
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could have a large impact on planetary evolution and structure ( [Morales et all 120091 ) 



For crystalline structures, the motion of the nuclei can be calculated in terms of 
oscillations about their equilibrium positions, i.e. phonon modes. In principle, phonon 
modes should be calculated self-consistently with electron excitations; in practice, for most 
conditions encountered in planets, the ion-thermal motion can be calculated from the 
phonon density of states at T = 0, □ and the electron-thermal contribution can often be 
ignored. f c comprises the ground state energy of the system as a function of compression, 
e c (p), and any configurational entropy that may be associa ted with different struct ural 



polymorphs or the entropy of mixing in an impure system ( IKittel fc Kroemer 



1980). One 



may calculate e c by setting up a configuration of nuclei, and finding the ground state of 
the Hamiltonian of the electrons. The configuration of the nuclei may be altered or relaxed 
under the net force they experience to find the lowest-energy structure, but they are often 
held fixed in likely structures. For well-defined crystalline structures, phonon modes can 
be predicted by calculating the electrostatic forces on the nuclei as each is displaced from 
equilibrium. Explicit calculation of the phonons is desirable below the Debye temperature, 
where the zero-point motion of the nuclei and the freezing out of modes may have a 



significant effect on the EOS ( Swift et a 



20011 ). 



Thermal excitation of the electrons can be calculated from the band structure, which 
is the set of eigenstates of the electronic Hamiltonian. For sufficiently high temperatures, 
the eigenstates must be calculated self-consistently with their population, but accurate 
calculations can be made to ~eV temperatures using the T = band structure. 

The key physics, and thus the source of limitations in the accuracy and validity of 
the EOS, is the electronic structure calculation. The challenge is in representing the 



2 A11 temperatures here are defined with respect to absolute zero. 
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multi-fermion nature of the electrons accurately, but with a method that is computationally 
tractable for real materials. Path integral methods, such as path-integral Monte-Ca rlo 



20091 ). 



(PIMC), can be used for direct, rigorous calculations of multi-particle states (IMilitzer 
However, these methods are practicable mostly for low-Z materials, and have not been 
applied systematically (or, generally, at all) to compositions of matter relevant to ice and 
rocky planets. There is, however, no reason to suspect that these techniques are unsuitable 
under exoplanet conditions per se. 



Although TFD-based electronic structure theory (ISalpeter fe Zapolskyill967D is often 



regarded as an adequate treatment for pressures over ~10 TPa (for instance 



2007 



Grasset et al 



Seager et al 



20091 ). it doe s not capture th e effects of electronic shell structure which 



are expected at high pressures ( jLiberman 



19791 ) . However, shell effects such as pressure 



ionization are captured by the electronic structure techniques considered here. 



Most theoretica 



1964 



Kohn fc Sham 



EOS c alculations us e dens ity functional theory (DFT) ( jHohenberg &: Kohn 



1965 



Perdew 



1992 



19941 ) and its variants to treat exchange and 



correlation between the electrons via functionals of the electron density, calibrated to 
reproduce the same system energy as techniques that treat the electron wavefunctions more 
directly. DFT functionals are calibrated against calculations of idealized electron gases, 
which may be performed up to arbitrarily high densities: they should be no less accurate 
for the relatively modest absolute compressions occurring in the cores of exoplanets than 



at p — 0. DFT calcu 
percent (for instance 



ations typically 



Swift et al 



2001 



r eproduce the o bserved mass density to within a few 



Ackland 



20021 ) . which is an observable discrepancy 



when compared with dynamic or quasistatic loading measurements. Calculations can be 
corrected by adjusting the energy to reproduce the mass density at zero pressure, which is 
most accurately measured. The resulting EOS - which we have termed 'ab fere initio' to 
distinguish them from uncorrected ab initio EOS - then typically reproduce high pressure 
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measurements to within their uncertainties ( ISwift et all 1200 ll . 120071 ) 



In many situations, electrons closer to the nuclei are not affected by changes in the 
compression or temperature, and the states of the other electrons may be calculated 
more efficiently by subsu ming the inner ele ctrons into a modified nuclear potential: the 



pseudopotential method (jPayne et al 



19921 ) . Pseudopotential calculations can become 
inaccurate when the nuclei are compressed sufficiently closely together. The validity can be 
checked by comparing against all-electro n calculations, and may be accurate to severalfold 



compressions or pressures of terapascals ([Morales et al 



2009|). 



With these caveats, pseudopotential techniques and the DFT construct in general are 
not inherently unsuitable for predicting EOS in the planetary structure regime. However, 
detailed calculations of many-species compounds are computationally expensive, and 
predictions of polymorphic structures are sensitive to relatively small inaccuracies in 
the computational methods. It is highly desirable to compare EOS predictions against 
experimental measurements. 

For the most part, we have used previously- developed EOS based on experimental data 
or validated electronic structure calculations up to the lim it of available data, and blending 



into TFD predictions at high compressions (IHolian 



1984 ). These EOS incorporate thermal 



effects, allowing us to investigate the sensitivity of planetary structures to temperature 
profile. For Fe and Fe-Ni, we have compared predictions from such EOS with 3D electronic 
structure calculations extended to higher pressures than have been reported previously. 
Previous theoretical studies include very careful treatments of the EOS of Fe up to pressur es 



and temperatures representative o 


'the 


Stixrude & Cohen 


1995 




Alfe et al 


2000 



Earth's core fin particular 



Sola et al 



20091 ; 



Wasserman et al 



Bclonosliko et al 



2008 



1996 



and 



references therein), and it is very desirable to perform equivalent studies to higher pressures 
and temperatures necessary for the study of exoplanets. For our present purposes, it is 
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most useful to compare EOS for different compositions of matter constructed according 
to consistent prescriptions. The use of wide-ranging EOS constructed using optimized 
algorithms, and wide-ranging sets of electronic structure calculations made using a 
consistent method that reproduces TFD at extreme compressions, allows us to avoid any 
reliance on extrapolating using ad hoc functional forms such as Vinet and Birch-Murnaghan 
for pressure-density relations, which can give unquantified uncertainties outside the range 
of the fitting data. 



2.2. Dynamic loading experiments 

The canonical experimental technique for studying the properties at matter at high 
pressure is shock loading, using a variety of methods to induce a shock wave. Shock 
measurements of EOS are often performed relative to a reference material, but an attraction 
of shock loading is that experiments can in principle be configured to yield absolute 
measurements, if the shock is induced by the impact of a projectile with a target of the 
same material. Indeed, the pressure standards in static compression apparatus such as 
diamond anvil cells are ultimately calibrated against absolute shock measurements. 

Although the time scales in dynamic loading experiments are typically nanoseconds 
to microseconds, typical equilibration times for electrons and atomic vibrations are much 
shorter, so inferred states used for testing and calibrating EOS are in thermodynamic 
equilibrium and thus equivalent to quasistatic compression measurements made in presses 
such as diamond anvil cells. Indeed, the difference in time scale between dynamic and 
quasistatic loading is less than the difference between the latter and planetary ages. 
Although EOS measurements are in thermodynamic equilibrium with respect to a given 
phase of matter, the time dependence of phase transitions must be considered: phase 
changes often occur with a significant degree of superpressurization or superheating 
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in dynamic loading experiments and indeed in quasistatic loading, compared with the 
equilibrium phase boundary. The effects of time dependence are also evident as hysteresis 
in the location of the phase change on loading compared with unloading. 

In shock loading, the entropy increases with compression, so the temperature rises 
faster with compression than it does along an isentrope. The Griineisen parameter and heat 
capacity of a material, which may be estimated theoretically or experimentally, can be used 
to predict the EOS away from the locus of shock states. Measurements along the principal 



shock Hugoniot hav e been used in t 



and the cold curve ( IBushman et al 



lis wa y for many years to estimate the principal isentrope 



19931 ). While a direct measurement of a relevant state 
is preferable, shock-derived EOS are likely to be adequate for exoplanets over a wide range 
of pressures. Furthermore, if a theoretical EOS is validated by shock experiments, this 
provides reasonable confidence that the EOS is valid at lower temperatures. A more serious 
limitation with shock experiments is that, for a given material and starting state, there is a 
limit to the compression that can be achieved by the passage of a single shock, and therefore 
a limit to the range of compressions that can be deduced along isentropes. However, this 
compression (discussed further below for different compositions of matter) is adequate to 
explore matter into the massive exoplanet regime. Shock measurements have been made 



to pr essures of up to 10 TPa using radiation from underground nuclear explosions (iRagan 
19841 ). 



The compression limit can be circumvented through the use of multiple shocks, which 
induce lower temperatures than a single shock to the same pressure. The ultimate limit 
for an infinite sequence of infinitesimal shocks is a ramp compression. Ramp compression 
of condensed matter was first demonstr ated using the expanding reaction products from 



a chemical explosive (iBarnes et al 



19741 1 . an d has subsequently been demonstrated using 



pulsed magnetic fields at the Z accelerator (IReisman et al 



20011 ) and high power lasers. 



13 



With lasers, the pulse energy may be used to vaporize a foil which then loa ds the sample 



20041 ) . the pulse may 



in an analogous way to the chemical reaction products (lEdwards et all . 
ablate the sample direc tly with an intensity history designed to induce a ramp in pressure 



( jSwift fc Johnson 



20051 ). or the sample may be driven from a hohlraum with a power 
history chosen to control the hohlraum temperature history such that a ramp is induced. 
Ramp loading using the last variant has been demonstrated to ~lTPa on the OMEGA 



laser ([Bradley et al 



20091 ). with near-term experiments on Fe planned to reach 2TPa at 



the National Ignition Facility (NIF). Experimental techniques have thus been available to 
develop and test EOS under conditions relevant to exoplanets, and planned developments 
should provide the first direct measurements of matter under exoplanet core conditions. 



3. Mass-radius relations 

If the rotation rate is not extreme, self-gravitating bodies are close to spherical. 
Spherical structures are certainly an adequate starting place for studies of exoplanets. 
The condition for isostatic equilibrium is that the stress induced by pressure variations is 
balanced by the gravitational acceleration g(r): 

gradp(r) = —p(r)g(r). (2) 

For Newtonian gravitation, by Gauss' theorem, g{r) can be expressed in terms of the mass 
enclosed within a given radius m(r): 

. , Gm(r) , . 

9(r) = —t 1 (3) 

so 

dgijr) _ c f 1 dm ( r ) 2 m ( r ) \ / 4 x 
dr \r 2 dr r 3 J 

m(r) can be calculated simply from the distribution of mass density, 

nr 

m(r) = 471 r' 2 p{r') dr' (5) 



o 
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or 

= 4nr 2 p(r). (6) 

dr 

When the total mass M, outer radius R, and surface composition are known with 
reasonable confidence, and the objective is to determine an aspect of the internal structure 



such as the core radius, the equations can be integrated from 
the known p(R) = (hence p(R)) and g(R) = GM/R 2 Jswift 



;he su rface inward, starting at 



20091 ) . The internal structure 



parameter is found as the solution of a shooting problem to ensure that m(r) — > as r — > 0. 
To calculate generic mass-radius relations for a uniform composition, it is most efficient to 
start at the center, where g(0) = and m(0) = 0, choose the core pressure p c = p(0), and 
integrate outward until p = 0. In this way, the mass-radius condition is obtained in terms 
of p c as an independent parameter: {R, M}(p c ). 

These solution methods are unwieldy when it is desired to constrain the overall 
composition as a mixture of components, such as a fixed ratio between Fe, rock, and H/He. 
In this case, the isostatic structure may be found more conveniently by starting with 
the desired masses of each component, and solving the system of equations as an elliptic 
problem. A numerical solution may be found iteratively by taking an approximate solution 
(which may initially be concentric shells of material at their p = density) and improving 
it by calculating the net stress gradient at each radius 

gradr = gradp + pg (7) 

and 

dp _ dpdp_ 
dr dp dr 

(ignoring the relatively slow change of g with deformation, which is valid near equilibrium), 
and hence estimating a change in radius for each shell such that the net pressure is zero. 
The radial variation was represented discretely by values at a finite series of radii. Multiple 
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iterations are necessary to allow the effect of a net pressure at some radius to be balanced 
through equilibration of the whole structure, but standard techniques may be used for 



19891 ). 



solving elliptic equations efficiently, such as successive over-relaxation (IPress et all 
A series of solutions describing a mass-radius relation can be found efficiently by working 
upward in mass, and increasing the density at each radius by a constant factor, then 
calculating the modified equilibrium structure. We verified that this method of solution 
gave the same result as integration from the center to the surface, by comparison with 
structures comprising pure Fe. 

To close the system of equations, a p(p) relation is needed for the material. The choice 
of p(p) relation implies a choice for the entropy or temperature profile of the body. In the 
present study we usually chose an isentrope passing through a reasonable surface state, such 
as STP for rock or metal, or a few atmospheres pressure and cryogenic temperatures for 
compositions that are gaseous at STP. Isentropes were calculated by numerical integration 
of the relation 

---p(p,e) (9) 



de 
dv 



20081 ). For a few sample 



using a procedure valid for EOS of arbitrary form (jSwift 
compositions, structures were calculated for p(p) along the cold curve. The precision of p(p) 
curves necessary to be useful for exoplanet structures is lower than for solar system planets, 
and the difference in structures between the STP isentrope and cold curve was found to be 
negligible over the range of pressures considered, as shown below. 

Mass-radius relations were calculated for a series of compositions of matter relevant 



to (exo)planets, using as a baseline EOS develop ed from and for shock wave applications, 

in| 



in particular the SESAME library (IHoliar 



19841 ). for which the constituent assumptions 



and calibrations are reasonably well documented. SESAME EOS are tabulated functions 
{p, e}(p, T). Most SESAME EOS are tabulated over a wide enough range of states 
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to be valid (at least in principle, to some finite precision) for massive exoplanets and 
star formation. They were constructed to incorporate experimental data, mainly shock 
measurements, where available. Experimental measurements not used in the construction of 
the EOS, including shock states, release from shock, and ramp compression, are frequently 
compared against the SESAME EOS, so there is a relatively rich literature on experimental 
validation. Many of the EOS were constructed using the same theoretical and modeling 
approaches, so conclusions drawn from the use of one EOS are more likely to apply to other 
materials. When these EOS are compared with subsequent shock or ramp measurements at 
higher pressures, it is quite usual to find significant differences, i.e. inaccuracies in the EOS, 
but they are usually at the level of a few tim es the experimental uncertainty (for instance 



Hicks et al 



2005 



Knudson fe Desjarlais 



20091 ). which would not have much effect on the 



mass-radius relation. 

For all undifferentiated compositions, the low-mass limit behaves as 

M = ^R 3 p (10) 

and deviates as matter within the planet is compressed gravitationally above its zero- 
pressure density po- The deviations become substantial as the central pressure approaches 
the bulk modulus. As the mass increases further, an increasing proportion is compressed to 
significantly higher density, until eventually any increase in mass is accommodated entirely 
by density increase, leading to a maximum in radius. 



3.1. Fe and Fe-Ni 



For Fe, which is the dominant component of interest at the highest pressures, we are 



developing 



(Swift 



)OS better optimized to exoplanet core states, to complement NIF experiments 



2010J). Thermodynamically complete versions of these EOS were not available for 
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this study: only the cold curves. The effect of temperature was assessed by comparing 
structures calculated with p(p) along an isentrope compared with the cold curve, for 
thermodynamically complete Fe EOS, and was found to be negligible, in agreement with 



previous work using a simplified thermal contribution ( jSeager et al 



20071 ). The mass-radius 



relation derived from our new Fe cold curve, which was calculated with DFT using the 
loca^ 



19911 ), 



density approximation and Troullier-Martins pseudopoten tials (ITroullier fe Martin 



agreed well with the relations from SESAME EOS 2140 flBarnes fe Rood 



2150 flKerlev 



19731 ) and 



19931 ) (Fig. [TJ). FesNi was considered as a representative Fe-Ni composition, 



calculated using t 



(ISwift et al 



2001 



re sa me DFT method as for previous studies of intermetallic compounds 
200 71). The p seudopotential for Ni was used previously for constructing 



EOS for NiAl ( Swift et al 



20071 ). Comparing p(p) relations, FeaNi was significantly softer 
than Fe or a segregated Fe-Ni mixture of the same composition, but the mass-radius 
relations did not differ significantly. The mass-radius relation was also calculated for 



the principal isent rope 



for rocky planets (ISwift 



rom t he linear Griineisen EOS found previously to perform well 



20091 ); it deviated from the other Fe relations for masses above 
approximately twice that of the Earth (Mg), demonstrating the inaccuracy of mass-radius 
relations deduced from EOS extrapolated from low-pressure measurements. Mass-radius 
relations deduced from the DFT cold curves were very similar to those from the SESAME 
EOS up to around 1 R/g or 3 Mg (~1 TPa), above which the properties calculated with DFT 
became progressively softer. This observation highlights the importance of ramp-loading 
experiments to several terapascals. With the exception of the linear Griineisen EOS, these 
relations are similar, but not identical, to relations used recently in the in terpretation o f 



exoplanet CoRoT-7b (jLeger et al 



2007 



Valencia et al 



20091 1 and for Fe mass-radius relations (IFortney et al 



20101 ) . Interestingly, the relation by Valencia et al follows the linear 



Griineisen curve up to around 2Mg, then approaches the T = DFT relations above 
10 Mg; the Fortney et al relation is stiffer (larger radius) up to around lMg, then follows 
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the relation for SESAME 2150 very closely. The mass-radius relations derived from the 
SESAME EOS were denser above 1M E . (Figs ED and H) 



For Fe, the maximum shock compression deduced from the SESAME EOS was a 



of aro und 5.2, calculated using a solution method valid for general forms of the EOS (ISwift 



actor 



20081 ). At the same mass density on the principal isentrope, the pressure is approximately 
llTPa. Thus single shock experiments, corrected appropriately for thermal effects, can 
plausibly explore planetary pressures to this regime. Impact-induced shock experiments, 
with the projectile driven by gas, chemical propella nt, or high exp losive, have been 



19661 ) . At high er pressures 



19841 ) 



performed on Fe to around 1 TPa by several groups ( Ivan Thiel 

data are very rare. Nuclear-driven shock impedance mismatch measurements ( iRagan 
reached a pressure which we estimate to be 4.57 ± 0.26 TPa. Further measurements above 
1 TPa are desirable. 



3.2. Rocks 

High quality theoretical and experimental studies have been performed for mantle 



Earth's core-man 


e boundary, and the effect of Fe ions 


Caracas & Cohen 


2008 




Stixn 


de k, Litheow-Bertelloni 


Stixrude & Litheow-Bertelloni 


2011 


). As compositions 



(see 


Drummond & Acklandl 


2002 


2010; 


Wentzcovitch et al 


2010 





in rocks, mass-radius relations were calculated for SiO?, M gO, and basalt. Si0 2 was 



represented by SESA ME EOS 7383 



SESAME EOS 7460 ( Barnes fc Lvon 



including quasiharmonic phonons (IDrummond Sz Ackland 



Joh nson fe Lyon 



19841 ) . MgO was represented by 



1988 a|) , and also by a DFT-derived two ph ase EOS 



2002 



Luo et a 



20041 ). but 



extending only to ~0.4TPa. The mass-radius relations for the different MgO models were 
indistinguishable, demonstrating insensitivity of the mass-radius relation to the inclusion of 
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Fig. 1. — Mass-radius relations calculated for different Fe EOS and for FesNi. 
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Fig. 2. — Predicted variation of core pressure with mass for different Fe EOS and for FesNi. 
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phase tran sitions, for current exop lanet applications. Basalt was represented by SESAME 



EOS 7530 (Barnes fc Lvon 



1988bl ). 'Basalt' refers to the composition (by number of atoms: 
60.11% O, 18.26% Si, 5.96% Al, 4.01% Ca, 3.48% Fe, 3.39% Mg, 2.18% H, 1.56% Na, 0.55% 
Ti, 0.38% K, 0.07% P, 0.06% Mn) and initial state: the EOS was derived from shock data 
including phase changes, and was blended into TFD predictions at high pressures. We also 
compared the linear Gruneisen E OS for basa lt found previously to give better than expected 



structures for the rocky planets ((Swift 



20091 ): its mass-radius relation was similar to that 
of the SESAME EOS to several tenths of a terapascal, and then diverged significantly. 
The mass-radius and pressure-radius relations were similar but distinguishable on scales 
extending to exoplanets, giving an indication of the possible uncertainty in inferences about 
rocky planet structure in the absence of constraints on the composition: the mass-radius 
relation does show sensitivity to the composition, but the sensitivity is small compared 
to current un c ertain ties in exoplanet observables. Recent mass-radius calculations by 



Valencia et all (120101 ) for Mg-silicate p 



relation; slightly older calculations by 



anets were a 



most i ndistinguishable from the MgO 



Fortney et all ( 120071 ) using a linear mixture of 38^ 



Si0 2 , 25% MgO, 25% FeS, and 12% FeO followed our calculation for basalt (SESAME) 
very closely to several tens of Mg. (Figs [3] and HI) 



3.3. Ices 



A b initio treatments of mixtures of molecules of low molecular weight, usually called 
ices (regardless of being in the solid or fluid state), are complicated and computationally 
expensive because of the possibility of segregation of species and the formation of a wide 
variety of chemical species. At issue here is the uncertainty of the equilibrium concentration 



of the different molecules. Such calculations 



numbers of thermodynamic states (jChau et al 



rave been performed for relatively small 



201 ll ) . and will in future allow wide-ranging 
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Fig. 3. — Mass-radius relations calculated for SiC>2, MgO, and basalt, in cluding a simplified, 
linear Griineisen EO S (LG), and comparing with previous calculations (IFortney et all 120071 ; 
Valencia et aJhoiff l. 




Fig. 4. — Predicted variation of core pressure with mass for SiC>2, MgO, and basalt, including 
a simplified, linear Griineisen EOS (LG). 
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EOS to be predicted as a function of composition. A complementary approach is to 
calculate the equilibrium chemical composition given the thermodynamic state and atomic 
composition, using measured or estimated thermodynamic potentials for each possible 
chemical specie. Chemical equilibrium calculations have been used to predict the EOS 
of the reaction products of chemical explosives up to temperatures of several hundred 
giga pascals and tempe r atures of ^eV, includi ng; many species in common with planetary 



ices (IBastea et al 



2007 



Mulford k Swift 



20061 ). However, the thermodynamic potentials are 



not necessarily well-constrained for the hot, compressed states occurring within exoplanets. 

As compositions representative of the variations in ice giants, mass-radius relations 
were calculated for H 2 0, NH 3 , and CH4, for isentropes chosen to pass through a condensed 
phase so that the planet surface could be defined as zero pressure. The calculated radius 
is thus a lower bound for planets with a thick atmosphere, as is likely to be the case for 
those orbiting close to their star. Several EOS were available for H?Q ; SESAME EOS 7154 



was the most appropriate for these pressures (IJohnson &: Lyon 



1990h. The STP isentrop e 



was used. The only SESAME EOS available for NH3 was 5520 (IHaar fc Gallagher 



1978|), 



fitted to National Bureau of Standards (NBS) gas phase data, and was tabulated to a 
maximum mass density of 0.765 g cm -3 , which is barely greater than typical for liquid 
ammonia. An empirical linear Griineisen EOS was constructed from shock wave data on 



20111 ). using an approximate 



liquid ammonia at an initial temperature of 203 K (jSwiftl 
treatment for off-Hugoniot states. The corresponding mass-radius calculation should be 
regarded with caution as it is not constrained by experimental data for pressures above a 
few tens of gigapascals. The mass-radius relation and central pressure were in reasonable 
agreement with predictions using SESAME EOS 5520 up to a mass ~0.1M^ above which 
the isentrope from the tab ular EOS bec ame unusable. For CH 4 , three SESAME EOS were 



available: 5500 and 5501 (iKerley 



19801 ) . both tabulated to 2.5 g cm and respectively 



with a Maxwell construction and van der Waals loops in the liquid-vapor region; and 
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5502 ( jJohnson 



1984j ). tabulated to 0.47 g cm -3 and based on NBS gas phase data. The 
mass-radius relation was calculated for 5500, which is the most relevant for planetary 
structures, for an isentrope passing through 70 K at zero pressure. The difference between 
H2O and CH4 was considerable: it is difficult to infer a composition for an ice giant from 
the mass-radius measurement alone to a factor ~2 in mass or ~25% in radius, without 



further constraints on the composition. (Figs and El} 



3.4. Hydrogen 



Jovian planets are thought to be composed predominantly of a mixture of H with 
approximately 5-10% He by number of atoms. 

H is theoretically interesting because quantum mechanical uncertainty in the position of 
the proton is more important than for other elements, and may affect the EOS. Significant 



work has been devo 



1992 



; ed to predicting the EOS of H and H/He mixtures ( jSaumon &: Chabrier 



Delaney et al 



2006 



Morales et al 



20091 ). These H/He mixtures are challenging 



because even the relatively small amount of helium involved is thought to greatly impact 
the EOS. S pecifically, the insu l ator-to-metal transit ion that occurs at o(100) GPa pressures 



in pure H ( [Morales et al 



2010 



Lorenzen et al 



20101 ) is thought to be suppres sed to even 



higher pressures and temperatures with the inclusion of trace amounts of He (jHensel 



1999|). 



Furthe rmore, helium is tho ught to segregate (i.e. form droplets) in the metallic region of the 



planet ([Morales et al 



20091 ). where the resulting He rain is thought to produce heat as the 
droplets fall down the gravitational well. This effect could have a disproportionately large 
effect on the mean EOS and on the differentiation of a planet by gravitational separation. 

We did not find a H/He EOS of similar pedigree to the other E OS above, so the 



mass-radius relation was calculated for pure H (SESAME EOS 5251: 



Albers fe Johnson 
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Fig. 5. — Mass-radius relations calculated for H2O, NH3, and CH4 ices. 
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mass (M Earth ) 

Fig. 6. — Predicted variation of core pressure with mass for H 2 0, NH 3 , and CH 4 ices. The 
crossing of the CH4 by the NH3 curve is at a mass density much greater than the shock data 
used to construct the NH 3 EOS, and should not be regarded as a meaningful prediction. 
The crossing of the H 2 and NH 3 curves is in a regime where the EOS should be adequately 
accurate. 
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( 119821 )) as a limiting case of low density. This EOS was isotopically scaled from S 



EOS 5263 for D 2 , which is wide-ranging and includes dissociation and ionization ( IKerley 



ISAM 



19721 ). and has been found to give good agreement with static and shock compression data. 
As with the ices, the isentrope was chosen such that a well-defined planetary surface was 
calculated, and is thus a lower bound on the radius of planets with a thick atmosphere. 
For the EOS used, the isentrope passed through 10 K at zero pressure, and exhibited a 
kink at around 100 caused by metallization at around 200 GP a. Interestingly, t he 



resulting mass-radius relation is quite close to the TFD prediction (ISeager et al 



3) 



2007). Our 



calculation has a maximum in radius, 79000 km at a mass of 11 x 10 16 kg. This result is 
slightly higher in mass and lower in radius than the TFD prediction. (Figs [7] and [HI) 



3.5. Fe/basalt 



Mass-radius relations were calculated for differentiated two-component planets 
comprising an Fe core and a basalt (composition) mantle, with fixed mass ratios, using the 
elliptic method of solution. Fe and basalt were represented by SESAME EOS 2150 and 
7530 respectively. Because of its low density, the addition of an outer basalt layer gave a 
disproportionate c hange in the radius, to a greater extent than indicated in previous work 



( ISeager et al 



20071 ). (Fig. El) 



4. Application to exoplanets 

For comparative purposes, we chose to represent each type of material by a single 
mass-radius curve corresponding to a single EOS: Fe-Ni by the electronic structure cold 
curve for FesNi (which is similar to the Fe SESAME 2150 curve but extends to higher 
pressure), rock by SESAME 5530 for basalt, ice by SESAME 7154 for H 2 0, and H/He by 
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Fig. 7. — Mass-radius relation calculated for H. 
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Fig. 8. — Predicted variation of core pressure with mass for H. 
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Fig. 9. — Mass-radius relations calculated for differentiated Fe/basalt planets, in terms of 
the ratio of basalt to Fe. 
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SESAME 5251 for H. Plotting the planets of the solar system, they fall in the expected 
places with respect to the curves: the rocky planets between Fe-Ni and rock with Mercury 
closest to Fe-Ni and Mars closest to rock, Jupiter and Saturn close to the H/He curve, and 
Uranus and Neptune close to the CH4 curve. Minor planets Pluto and Eris are also shown, 
lying on the icy side of the basalt cur ve. We also compar e exoplanets at extremes of siz e 



and mass: 'sup er-Earths' Kepler-lOb (|Batalha et alll201ll ). CoRoT-7b ( 



and GJ 1214b (ICharbonneau et al 



20081 ) and CoRoT-3b flDeleuil et al 



Leger et al 



2009) and the 'super- Jupiters' HAT-P-2b flLoeillet et al 



2009), 



20081 ). whose core pressures are representative of 



states it would be useful to explore in future experiments on material propert ies, and 



the an omalously large CoRoT-2b. (Figs [10] and [TTJ exoplanet parameters from 

fboiit u 



Schneider 



Measurements of exoplanet Kepler-lOb constrain its radius to 4.56l}'2gM£ and 
mass to 1.416^o'g3g Re- The nominal value lies on our mass-radius relation for Fe-basalt 
differentiated planets with 2/3 of the total mass in the core. With this structure, the 
central pressure would be 2. 5^93 TPa. Considering two-layer differentiated structures at 
the la level, the planet could have a core mass fraction between 25% and 95%, with central 
pressure between about 1.5 and 3.7 TPa. 

Measurements of exoplanet C0R0 T- 7b constrain its radius to 1.68 ± 0.09 Rg and 



mass to 4.8 ± 0.8 Mg (ILeger et al 



20091 ) . though other resear chers have deduced different. 



and inconsistent, masses from the same Doppler shift data (IHatzes et al 



2010 



Pont et al 



2011 



Ferraz-Mello et al 



20111 ) . Given the current discrepancy in mass deduced by different 
groups, the validity of deduced compositions and central pressures is unclear, but it is 
informative to assess the values and un certainties derived from the quoted uncertainties in 



mass and radius. The nominal value by ILeger et all ( 120091 ) lies very close to our mass-radius 
relation for basalt, for which the central pressure would be 0.751^2 TPa. Considering 
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Fig. 10. — Mass-radius relations for different classes of material, compared with planets, 
moons, and exoplanets (red crosses). 
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Fig. 11. — Mass-radius relations for different classes of material, compared with selected 
exoplanets. 
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two-layer differentiated structures at the la level, the planet could be Fe/rock with a 
core of radius up to 4500 km, or rock/ice with a surface layer of (for example) H 2 up 
to 2500 km thick. The central pressure is calculated to increase rapidly with Fe, giving 
0.8loiTPa. For three-layer differentiated structures, the ice layer would be thicker as 
the core size increased. The uncertainty in deduced composition, core radius, and central 
pressure is nonlinear with respect to the uncertainty in mass and radius. Measurements of 
Doppler shift in the stellar spectrum and dimming in the stellar brightness are essentially 
independent, so we calculate uncertainty contours in radius and central pressure and 
smooth them to infer the la uncertainty. The planet is v ery unlikely to be Fe-only, in 



agreement with the previous conclusions ( iLeger et al 



20091 ). In fact, CoRoT-7b is likely to 



have a smaller core in proportion to its size than do the terrestrial planets. The presence 
of ice is thought to be unlikely because of the proximity of the planet to its star, CoRoT-7, 
resulting in high surface temperatures. The nominal mass-radius measurement is consistent 
with a composition of rock only; given the uncertainties in the mass-radius measurement 
and relations, a metallic core cannot be ruled out. 



Recent measu rements of exoplanet G 



1214b give a radius of 2.65 ± 0. 



1 Rp and mass 



20111), lying 



of 6.55 ± 0.9 M E flCharbonneau et all 120091 ) or 6.45 ± 0.9 M E (ICarter et all 
between the mass-radius curves for H 2 and CH 4 . The closest analogue in the solar system 
is the icy giants; the mass and radius of GJ 1214b are consistent with a higher proportion of 
H2O, or with a proportionately larger amount of rock or Fe, than Neptune or Uranus. The 
observations do not of themselves suggest an 'ocean planet' as the planet's density see ms 



implausibly low. Furthermore, as has been pointed out previously (lAdams et al 



20081 ). a 



mean density consistent with a 'water world' could arise from a rocky interior with an H/He 
atmosphere. 



Recent parameters for the super-massive exoplanet HAT-P-2b are a mass of 
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9.09 ± 0.24Mj and a radius of 1.157^;^ Rj flPal et aJl2009t ). These parameters are 
consistent with a massive hydrogen- rich planet, as concluded previously. Proportionately, 
the radius is smaller with respect to our predicted mass-radius curve for hydrogen than are 
Jupiter or Saturn, suggesting a larger proportion of rock or Fe. A pure hydrogen planet of 
the nominal mass would have a radius of 1.12 MJ and a central pressure of 214 TPa. The 
actual proportion of rock that would be consistent with the observed mass and radius is 
difficult to estimate using the available EOS, because they roll over at lower masses; future 
studies will include theoretical prediction and experimental investigation of the EOS under 
relevant conditions. 

The super-ma ssive exoplanet CoRoT-3b has a mass of 21.7 ± 1 Mj and 1.01 ± 0.07 Rj 



(IDeleuil et al 



20081 ) . lying close to our mass-radius curve for H. 



Our mass-radi us relation for H co nfirms the anomalously large radius of CoRoT-2b. As 



proposed recently (ILeconte et al 



20091 ). a plausible explanation is that this planet is young 
and hot, from accretion, tidal heating, or internal processes. The radius would be larger for 
a warmer isentrope, but is relatively insensitive to temperature. Another potential source 
of uncertainty, particularly for warm bodies with a significant proportion of volatiles, is 
the sensitivity to definition of the outer surface of the body. If the low-density gas in the 
outer region of the atmosphere has a higher opacity than expected, the planet would appear 
larger. 

Using the location of a given planet with respect to the mass-radius relation for 
different compositions, it is possible to estimate the pressure in the center of the core from 
the corresponding mass-pressure relation (Fig. fl~2l) . Thus the central pressure of Kepler-lOb 
is likely to be between 1.5 and 2.7 TPa, and that of CoRoT-7b is likely to be between 0.7 
and 2 TPa, which are accessible using planar ramp-loading experiments at the OMEGA 
and NIF facilities. The core pressure in GJ 1214b is approximately lTPa: the same mass 
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of H2O gives a pressure of around 1.5 TPa, and that mass of CH4, 0.5 TPa. It would be 
most useful to improve EOS of ice compositions into this regime, by QMD simulations and 
dynamic loading experiments; again, planar ramp experiments at NIF should be able to 
access relevant states. The core pressure in HAT- P- 7b is approximately 20 TPa, which may 
be accessible with subsequent NIF experiments using convergent compression. The core 
of CoRoT-3b is at the highest pressure known for non-stellar matter: 1900 TPa according 
to the EOS used. Such pressures are achieved in inertial confinement fusion experiments, 
but will require very significant experimental developments before EOS experiments can be 
performed. 



5. Conclusions 

Mass-radius relations were derived using EOS widely used for, and thus validated 
against, dynamic compression experiments. The results were similar to relations used 
previously to infer planetary structures, but were not identical. 

Mass-radius relations for Fe were inferred using EOS calculated to higher pressures 
than reported previously, using density functional theory. The difference between the cold 
curve and STP isentrope had a negligible effect on the mass-radius relation in the exoplanet 
regime. The effect of Ni in Fe was significant for the EOS, but negligible for the mass-radius 
relation of exoplanets. In a differentiated planet, the radius was more sensitive to small 
proportions of rock than has been concluded previously. 

These results support the conclusion that CoRoT-7b, with measurements analyzed by 



Leger et all (120091 ) . is 'Earth-like' in that it is likely to have a rocky layer over an Fe-Ni core, 
though the present uncertainty in the exoplanet's mass is too large to constrain the mass 
of Fe significantly: there may be very little. GJ 1214b lies between the curves for H 2 and 



-38 - 




0.001 0.01 0.1 1 10 100 1000 10000 

mass (M Earth ) 

Fig. 12. — Mass-pressure relations for different classes of material. 
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CH 4 , which are compositions where improvements to the EOS are particularly desirable. 
The composition cannot be constrained meaningfully from present data, and suggestions 
that it might be an 'ocean planet' are speculative. The present results also support the 
conclusion that HAT-P-2b is a gas giant, though the density is proportionately larger than 
the gas giants in the solar system, indicating a higher proportion of matter of higher atomic 
number. The mass-radius relations imply core pressures of 0.75 ± 0.15TPa for CoRoT-7b, 
a range range 0.5-1.5 TPa for GJ 1214b, over 210 TPa in HAT-P-2b, and around 1900 TPa 
in CoRoT-3b. 

Based on the constituent physics and approximations, we conclude that current 
theoretical techniques for electronic structure, with care in the use of pseudopotentials, 
should be as suitable for predicting EOS in the exoplanet regime as they are at low 
pressures. Shock measurements are relevant in terms of the compressions that can be 
explored, though interpretation is involved to account for the elevated temperatures. Ramp 
loading techniques are very appropriate, and planned NIF experiments should provide data 
on Fe in the region where DFT mass-radius relations deviate from SESAME EOS. In both 
experiment and theory, the limiting cases of Fe (or Fe-Ni) and H are the subject of active 
research. Rocks and ices are more complicated in terms of composition and structure, 
and there is clearly a need for future research into the properties of these materials at 
compressions relevant for planets and exoplanets. 
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